home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT3.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  11KB  |  386 lines
  1. //$$ newmat3.cxx        Matrix get and restore rows and columns
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5.  
  6. #include "include.hxx"
  7.  
  8. #include "boolean.hxx"
  9.  
  10. typedef double real;                 // all floating point double
  11.  
  12. #include "newmat.hxx"
  13.  
  14. #include "newmatrc.hxx"
  15.  
  16. //#define REPORT { static ExeCounter ExeCount(__LINE__,3); ExeCount++; }
  17.  
  18. #define REPORT {}
  19.  
  20. //#define MONITOR(what,storage,store) \
  21. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  22.  
  23. #define MONITOR(what,store,storage) {}
  24.  
  25.  
  26. // Control bits codes for GetRow, GetCol, RestoreRow, RestoreCol
  27. //
  28. // LoadOnEntry:
  29. //    Load data into MatrixRow or Col dummy array under GetRow or GetCol
  30. // StoreOnExit:
  31. //    Restore data to original matrix under RestoreRow or RestoreCol
  32. // IsACopy:
  33. //    Set by GetRow/Col: MatrixRow or Col array is a copy
  34. // DirectPart:
  35. //    Load or restore only part directly stored; must be set with StoreOnExit
  36. //    Still have decide  how to handle this with symmetric
  37. // StoreHere:
  38. //    used in columns only - store data at supplied storage address, adjusted
  39. //    for skip; used for GetCol, NextCol & RestoreCol. No need to fill out
  40. //    zeros.
  41.  
  42.  
  43. // These will work as a default
  44. // but need to override NextRow for efficiency
  45.  
  46.  
  47. void GeneralMatrix::NextRow(MatrixRowCol& mrc)
  48. {
  49.    REPORT
  50.    if (mrc.cw*StoreOnExit) { REPORT this->RestoreRow(mrc); }
  51.    if (mrc.cw*IsACopy)
  52.    {
  53.       REPORT
  54.       real* s = mrc.store + mrc.skip;
  55.       MONITOR("Free   (NextRow)",mrc.storage,s) 
  56.       delete [mrc.storage] s;
  57.    }
  58.    mrc.rowcol++;
  59.    if (mrc.rowcol<nrows) { REPORT this->GetRow(mrc); }
  60.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  61. }
  62.  
  63. void GeneralMatrix::NextCol(MatrixRowCol& mrc)
  64. {
  65.    REPORT                                                // 423
  66.    if (mrc.cw*StoreOnExit) { REPORT this->RestoreCol(mrc); }
  67.    if (mrc.cw*IsACopy && (mrc.cw*StoreHere==0))
  68.    {
  69.       REPORT                                             // not accessed
  70.       real* s = mrc.store + mrc.skip;
  71.       MONITOR("Free   (NextCol)",mrc.storage,s) 
  72.       delete [mrc.storage] s;
  73.    }
  74.    mrc.rowcol++;
  75.    if (mrc.rowcol<ncols) { REPORT this->GetCol(mrc); }
  76.    else { REPORT mrc.cw -= (StoreOnExit+IsACopy); }
  77. }
  78.  
  79.  
  80. // routines for matrix
  81.  
  82. void Matrix::GetRow(MatrixRowCol& mrc)
  83. {
  84.    REPORT
  85.    mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=ncols;
  86.    mrc.store=store+mrc.rowcol*ncols;
  87. }
  88.  
  89.  
  90. void Matrix::GetCol(MatrixRowCol& mrc)
  91. {
  92.    REPORT
  93.    mrc.skip=0; mrc.storage=nrows;
  94.    if (ncols==1 && (mrc.cw*StoreHere)==0)
  95.       { REPORT mrc.cw-=IsACopy; mrc.store=store; }           // not accessed
  96.    else
  97.    {
  98.       mrc.cw+=IsACopy; real* ColCopy;
  99.       if ((mrc.cw*StoreHere)==0)
  100.       {
  101.          REPORT
  102.          ColCopy = new real [nrows]; MatrixErrorNoSpace(ColCopy);
  103.          MONITOR("Make (MatGetCol)",nrows,ColCopy)
  104.          mrc.store = ColCopy;
  105.       }
  106.       else { REPORT ColCopy = mrc.store; }
  107.       if (mrc.cw*LoadOnEntry)
  108.       {
  109.          REPORT
  110.          real* Mstore = store+mrc.rowcol; int i=nrows;
  111.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  112.       }
  113.    }
  114. }
  115.  
  116. void Matrix::RestoreCol(MatrixRowCol& mrc)
  117. {
  118. //  if (mrc.cw*StoreOnExit)
  119.    REPORT                                   // 429
  120.    if (mrc.cw*IsACopy)
  121.    {
  122.       REPORT                                // 426
  123.       real* Mstore = store+mrc.rowcol; int i=nrows; real* Cstore = mrc.store;
  124.       while (i--) { *Mstore = *Cstore++; Mstore+=ncols; }
  125.    }
  126. }
  127.  
  128. void Matrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrMat(); }  // 1808
  129.  
  130. void Matrix::NextCol(MatrixRowCol& mrc)
  131. {
  132.    REPORT                                        // 632
  133.    if (mrc.cw*StoreOnExit) { REPORT RestoreCol(mrc); }
  134.    mrc.rowcol++;
  135.    if (mrc.rowcol<ncols)
  136.    {
  137.       if (mrc.cw*LoadOnEntry)
  138.       {
  139.          REPORT
  140.          real* ColCopy = mrc.store;
  141.          real* Mstore = store+mrc.rowcol; int i=nrows;
  142.          while (i--) { *ColCopy++ = *Mstore; Mstore+=ncols; }
  143.       }
  144.    }
  145.    else { REPORT mrc.cw -= StoreOnExit; }
  146. }
  147.  
  148. // routines for diagonal matrix
  149.  
  150. void DiagonalMatrix::GetRow(MatrixRowCol& mrc)
  151. {
  152.    REPORT
  153.    mrc.skip=mrc.rowcol; mrc.cw-=IsACopy; mrc.storage=1; mrc.store=store;
  154. }
  155.  
  156. void DiagonalMatrix::GetCol(MatrixRowCol& mrc)
  157. {
  158.    REPORT 
  159.    mrc.skip=mrc.rowcol; mrc.storage=1;
  160.    if (mrc.cw*StoreHere) 
  161.       { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); mrc.cw+=IsACopy; }
  162.    else { REPORT mrc.store = store; mrc.cw-=IsACopy; }     // not accessed
  163. }
  164.  
  165. void DiagonalMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrDiag(); }
  166.                               // 800
  167. void DiagonalMatrix::NextCol(MatrixRowCol& mrc)
  168. {
  169.    REPORT
  170.    if (mrc.cw*StoreHere)
  171.    {
  172.       if (mrc.cw*StoreOnExit) 
  173.          { REPORT *(store+mrc.rowcol)=*(mrc.store+mrc.rowcol); }
  174.       mrc.IncrDiag();
  175.       if (mrc.cw*LoadOnEntry && mrc.rowcol < ncols)
  176.          { REPORT *(mrc.store+mrc.rowcol)=*(store+mrc.rowcol); }
  177.    }
  178.    else { REPORT mrc.IncrDiag(); }                     // not accessed
  179. }
  180.  
  181. // routines for upper triangular matrix
  182.  
  183. void UpperTriangularMatrix::GetRow(MatrixRowCol& mrc)
  184. {
  185.    REPORT
  186.    int row = mrc.rowcol; mrc.skip=row; mrc.cw-=IsACopy;
  187.    mrc.storage=ncols-row; mrc.store=store+(row*(2*ncols-row-1))/2;
  188. }
  189.  
  190.  
  191. void UpperTriangularMatrix::GetCol(MatrixRowCol& mrc)
  192. {
  193.    REPORT
  194.    mrc.skip=0; mrc.cw+=IsACopy; int i=mrc.rowcol+1; mrc.storage=i;
  195.    real* ColCopy;
  196.    if ((mrc.cw*StoreHere)==0)
  197.    {
  198.       REPORT                                              // not accessed
  199.       ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
  200.       MONITOR("Make (UT GetCol)",i,ColCopy) 
  201.       mrc.store = ColCopy;
  202.    }
  203.    else { REPORT ColCopy = mrc.store; }
  204.    if (mrc.cw*LoadOnEntry)
  205.    {
  206.       REPORT
  207.       real* Mstore = store+mrc.rowcol; int j = ncols;
  208.       while (i--) { *ColCopy++ = *Mstore; Mstore += --j; }
  209.    }
  210. }
  211.  
  212. void UpperTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  213. {
  214. //  if (mrc.cw*StoreOnExit)
  215.   {
  216.      REPORT
  217.      real* Mstore = store+mrc.rowcol; int i=mrc.rowcol+1; int j = ncols;
  218.      real* Cstore = mrc.store;
  219.      while (i--) { *Mstore = *Cstore++; Mstore += --j; }
  220.   }
  221. }
  222.  
  223. void UpperTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrUT(); }
  224.                               // 722
  225.  
  226. // routines for lower triangular matrix
  227.  
  228. void LowerTriangularMatrix::GetRow(MatrixRowCol& mrc)
  229. {
  230.    REPORT
  231.    int row=mrc.rowcol; mrc.skip=0; mrc.cw-=IsACopy; mrc.storage=row+1;
  232.    mrc.store=store+(row*(row+1))/2;
  233. }
  234.  
  235. void LowerTriangularMatrix::GetCol(MatrixRowCol& mrc)
  236. {
  237.    REPORT
  238.    int col=mrc.rowcol; mrc.skip=col; mrc.cw+=IsACopy;
  239.    int i=nrows-col; mrc.storage=i; real* ColCopy;
  240.    if ((mrc.cw*StoreHere)==0)
  241.    {
  242.       REPORT                                            // not accessed
  243.       ColCopy = new real [i]; MatrixErrorNoSpace(ColCopy);
  244.       MONITOR("Make (LT GetCol)",i,ColCopy) 
  245.       mrc.store = ColCopy-col;
  246.    }
  247.    else { REPORT ColCopy = mrc.store+col; }
  248.    if (mrc.cw*LoadOnEntry)
  249.    {
  250.       REPORT
  251.       real* Mstore = store+(col*(col+3))/2;
  252.       while (i--) { *ColCopy++ = *Mstore; Mstore += ++col; }
  253.    }
  254. }
  255.  
  256. void LowerTriangularMatrix::RestoreCol(MatrixRowCol& mrc)
  257. {
  258. //  if (mrc.cw*StoreOnExit)
  259.    {
  260.       REPORT
  261.       int col=mrc.rowcol; real* Cstore = mrc.store+col;
  262.       real* Mstore = store+(col*(col+3))/2; int i=nrows-col;
  263.       while (i--) { *Mstore = *Cstore++; Mstore += ++col; }
  264.    }
  265. }
  266.  
  267. void LowerTriangularMatrix::NextRow(MatrixRowCol& mrc) { REPORT mrc.IncrLT(); }
  268.                                      //712
  269. // routines for symmetric matrix
  270.  
  271. void SymmetricMatrix::GetRow(MatrixRowCol& mrc)
  272. {
  273.    REPORT                                                //571
  274.    mrc.skip=0; int row=mrc.rowcol;
  275.    if (mrc.cw*DirectPart)
  276.    {
  277.       REPORT
  278.       mrc.cw-=IsACopy; mrc.storage=row+1; mrc.store=store+(row*(row+1))/2;
  279.    }
  280.    else
  281.    {
  282.       mrc.cw+=IsACopy; mrc.storage=ncols;
  283.       real* RowCopy = new real [ncols]; MatrixErrorNoSpace(RowCopy);
  284.       MONITOR("Make (SymGetRow)",ncols,RowCopy) 
  285.       mrc.store = RowCopy;
  286.       if (mrc.cw*LoadOnEntry)
  287.       {
  288.      REPORT